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Chapter  1 


INTRODUCTION 
1.1  Motivation 

New-generation  military  aircraft  are  being  developed  to  take  off  and  land  vertically,  resulting 
in  large  thermal  loads  on  the  pavement  surface.  These  Fast  transient  thermal  loads  produce 
a  rapidly  varying  temperature  profile  both  radially  and  through  the  slab.  These  repeated 
thermal  loads  can  lead  to  premature  deterioration  of  the  airfield  pavement  structure.  Ju  and 
Zhang1,2  give  a  detailed  account  on  this  issue.  Traditional  paving  materials  such  as  concrete 
and  asphalt  concrete  will  not  have  the  same  longevity  under  this  condition  of  repeated 
thermal.1-6  Accurately  predicting  this  transient  high-temperature  profile  is  crucial  and  a 
prerequisite  to  further  determining  the  thermal  stress  fields  in  the  material  design  of  this 
new  type  of  airfield  pavement  application. 

Different  approaches  can  be  applied  to  predict  temperature  fields  in  multi-layered  pave¬ 
ment  systems  under  climatic  conditions,  such  as  statistics-based  models,  numerical  ap¬ 
proaches,  and  analytical  methods.  Wang,  Roesler,  and  Guo'  present  an  overview  of  these 
various  approaches.  To  estimate  rapidly  changing  temperature  profiles  in  concrete  pavements 
subjected  to  fast  transient  thermal  loads,  numerical  or  analytical  approach  is  appropriate.  A 
numerical  approach — specifically,  an  explicit  finite  difference  method — was  employed  to  pre¬ 
dict  two-dimensional  (2-D)  axisymmetric  transient  high  temperature  field  in  Ju  and  Zhang.2 
The  main  advantage  of  this  method  is  that  it  can  easily  handle  temperature-dependent  ther¬ 
mal  properties  of  concrete,  such  as  density,  specific  heat,  and  thermal  conductivity.  However, 
an  extremely  small  temporal  step  size,  which  is  highly  dependent  on  the  spatial  step-size, 
must  be  chosen  to  ensure  computational  stabilities  and  hence  this  approach  is  generally  more 
time-consuming.  In  view  of  the  facts  that  the  thermal  properties  of  concrete  change  slowly 
when  the  temperature  increases,2  and  the  huge  amounts  of  heat  exhausted  from  aircraft 
engines  is  the  dominant  driving  force  for  this  problem,  a  rapid  analytical  solution  can  be 
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pursued  by  assuming  that  thermal  properties  of  each  layer  of  material  in  pavements  are 
constant. 

The  primary  objective  of  this  study  is  to  develop  easily  implemented  analytical  solutions 
for  predicting  the  rapidly  varying  pavement  temperature  profile  under  fast  transient  thermal 
loads.  To  meet  different  needs  for  acquiring  temperature  profile  information,  one-dimensional 
(1-D)  and  2-D  axisymmetric  temperature  fields  will  be  considered.  Time-dependent  temper¬ 
ature  profiles  change  only  with  depth  in  the  1-D  case,  and  the  profiles  change  both  vertically 
and  radially(horizontally)  in  the  2-D  axisymmetric  case.  The  main  advantage  of  these  an¬ 
alytical  solutions  is  that  they  lay  the  foundation  for  further  investigating  the  1-D  and  2-D 
axisymmetric  thermal  stress  fields  in  concrete  pavements  based  on  analytical  approaches. 

1.2  Overview 

The  overview  of  this  report  is  as  follows: 

In  Chapter  2,  1-D  temperature  fields  in  homogeneous  half-space  subjected  to  fast,  tran¬ 
sient,  thermal  loadings  is  investigated  first.  The  general  closed- form  solution  for  this  initial¬ 
boundary  problem  is  identified  through  literature  review.  Efficient  Gaussian-type  quadrature 
formulas  developed  by  Steen  et  al.  are  tested  and  recommended  to  numerically  resolve  the 
general  solution.  This  is  followed  by  the  study  of  1-D  temperature  fields  in  two-layered 
pavement  systems  subjected  to  high-temperature  transient  thermal  loadings.  Two  types  of 
solutions  are  derived  based  on  two  different  boundary  conditions,  namely,  specified  pavement 
surface  temperature  history,  and  specified  heat  flux  intensity  from  aircraft  engines,  respec¬ 
tively.  The  main  mathematical  tools  employed  in  deriving  analytical  solutions  in  these  cases 
are  Laplace  integral  transforms  (LT)  and  numerical  Laplace  inversion.  Some  model  calcula¬ 
tions  are  performed  to  demonstrate  the  derived  analytical  solutions.  This  chapter  concludes 
with  some  sensitivity  studies  that  investigate  effects  of  material  thermal  properties  and 
thickness  of  the  first  layer  on  temperature  profile  in  a  two-layered  pavement  system. 

In  Chapter  3,  a  2-D  axisymmetric  temperature  held  in  homogeneous  half-space  with  spec- 
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ified  surface  temperatures  clue  to  transient  thermal  loads  is  studied.  Analytical  solutions  are 
derived  respectively  using  two  different  methods:  one  based  on  separation  of  variables  (SV) 
and  HT  while  the  other  is  based  on  LT  and  HT.  Numerical  results  are  obtained  using  the 
derived  solutions  and  a  model  pavement  surface  temperature  history  extracted  from  Dr. 
Zhang’s  doctoral  dissertation,8  suggesting  that  a  combined  result  based  on  these  two  differ¬ 
ent  analytical  solutions  will  give  a  reasonable  approximation  of  the  pavement  temperature 
profile. 

In  Chapter  4,  a  summary  of  this  technical  report  is  presented. 
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Chapter  2 


1-D  TEMPERATURE  FIELD  IN  PAVEMENTS 

In  this  chapter,  1-D  rapidly  varying  temperature  profiles  in  pavements  subjected  to  fast, 
transient  thermal  loadings  is  investigated.  In  Section  2.1,  an  analytical  solution  for  such 
temperature  profiles  in  a  homogeneous  half-space  is  presented.  Analytical  solutions  for 
temperature  profile  in  two-layered  pavement  systems  are  systematically  studied  using  LT 
and  numerical  Laplace  inversions  in  Section  2.2. 

2.1  1-D  Temperature  Field  in  Homogeneous  Half-Space  Subjected  to  Fast, 
Transient  Thermal  Loadings 

2.1.1  Specified  Pavement  Surface  Temperature 

The  governing  equation  for  this  heat  conduction  problem  without  an  internal  heat  source/sink 
is  the  classic  1-D  heat  equation 

d  T  d2T 

— —  =  a——  for  0  <  z  <  oo  and  t  >  0  (2.1) 

at  oz 2 

where  a  =  thermal  cliffusivity  of  material  (m2/h). 

One  way  to  consider  rapidly  transient  thermal  loadings,  e.g.,  energy  projected  from 
vertical-take-off/landing  aircraft  with  fast  heating  rate  (say,  500  °C/min.  as  used  in  Ju  and 
Zhang1),  is  to  use  measured  transient  surface  temperatures  F(t)  (if  available)  in  the  area 
where  the  temperature  is  the  highest.  Mathematically,  the  following  initial  boundary  value 
problem  needs  to  be  solved 


dT .  . 

d2T 
adz 2' 

(■ z,t ) 

0  <  z  <  oo  and  0  <  t  <  oo 

T(0,t) 

=  m 

for 

z  =  0 

(2.2) 

T(z,  0) 

=  G(z) 

for 

t  =  0 
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The  analytical  solution  for  the  above  initial  boundary  value  problem  (Equation  2.2)  can 
be  obtained  using  the  method  of  odd  extension  discussed  in  Section  3.1  in  Strauss,9  or  by 
summing  up  solutions  of  two  simpler  initial-boundary-value  problems  outlined  on  page  64 
in  Carslaw  and  Jaeger.10  The  complete  solution  of  Equation  2.2  is 


T{z,t) 


(z+y)2 

4  at 


G{y)dy 


e  y2  dy 


provided  the  improper  integrals  in  Equation  (2.3)  converge. 

In  Ju  and  Zhang,1  F(t)  and  G(z)  take  the  following  forms 


m  =  T,(t) 

=  285  + 49.5  ln(i  + 0.00554) 
G(z)  =  T„ 

=  25 


(2.3) 


(2.4) 


(2.5) 


where  t  in  Equation  (2.4)  is  measured  in  seconds,  and  Tq  in  Celsius  degrees  in  Equation 
(2.5). 

Substituting  Eqs.  (2.4)  and  (2.5)  into  Equation  (2.3)  gives 

T{z' t]  =  h  L_ z  (‘  ■  ssy )  e^dy +-^T°  Sr  e^dy  (2  6) 

\/4  at 

which  is  in  agreement  with  Equation  (5)  in  Ju  and  Zhang.1 

It  is  noted  that  for  arbitrary  z  >  0  and  t  >  0,  the  improper  integral  in  Equation  (2.6) 
can  be  shown  to  converge  to  a  finite  value  by  using  the  Lebesgue’s  dominated  convergence 
theorem  from  real  analysis.  The  time-dependent  surface  temperature  described  by  Equation 
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(2.4)  is  plotted  in  Figure  2.1. 


Figure  2.1:  Time-dependent  Surface  Temperature  Due  to  Transient  High- 
temperature  Loadings 


Due  to  the  complexities  of  integrands  in  the  integrals  in  Equation  (2.6),  the  closed-form 
solution  of  integrals  in  Equation  (2.6)  is  hard  to  derive.  Thus  a  numerical  approximation  to 
Equation  (2.6)  is  employed  in  this  study. 

Steen  et  al.11  developed  efficient  Gaussian- type  integration  formulas  to  approximate 
integrals  of  the  forms 

r°°  2  r1  2 

/  e  x  f(x)dx  and  /  e  x  f(x)dx , 

Jo  Jo 
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and  their  formulas  are  given  as  follows 


1 


e  X'f(x)dx 


N 

2>/(**) 

k= 1 


M 

^2wkf(xk) 

k=  1 


(2.7) 

(2.8) 


where  N  e  {k  is  a  positive  integer  :  2  <  k  <  15};  M  e  {/c  is  a  positive  integer  :  2  <  k  < 
10};  weights  wk  and  abscissae  xk,  k  —  1,  2, . . , ,  N  or  M  are  listed  in  Steen  et  al.11 

To  apply  Equations  (2.7)  and  (2.8),  the  integrals  in  Equation  (2.6)  have  to  be  transformed 
into  the  standard  integral  form  J0°°  e~x~  f(x)dx  or  e~x"  f{x)dx.  This  can  be  easily  achieved 
using  a  change  of  variables  as  follows: 

Let  rj  =  -j=  and  y  =  rj  +  £;  then  the  improper  integral  in  Equation  (2.6)  becomes 


Ts  t 


4  ay2 


e~y2 dy  =  /  ( t 


4  a(v  +  O2 


d{ 


(2.9) 


on  the  other  hand,  if  we  let  y  =  r)£,  the  definite  integral  in  Equation  (2.6)  becomes 


q  f  e  ^e1-1  r'2^2  d£ 


(2.10) 


For  fixed  z  >  0,  t  >  0,  the  temperature  T(z,  t )  can  be  approximated  by  applying  Steen  et 
al.’s11  integral  formulas  to  Equation  (2.6).  To  investigate  the  effect  of  the  concrete  diffusivity 
coefficient  a  on  the  temperature  profile,  a  =  1.3  mm2/s  and  a  =  1.0  mm2/s  used  by  Ju 
and  Zhang1  are  adopted  in  this  study,  and  IV  =  15,  M  =  10  are  employed  in  Steen  et  al’s11 
integral  formulas.  For  the  sake  of  completeness  Tables  2.1  and  2.2  list  the  weights  wk  and 
abscissae  xk  used  in  Eqs.  (2.7)  and  (2.8),  respectively. 

Figure  2.2  plots  effects  of  concrete  diffusivity  coefficient  a  on  concrete  pavement  temper¬ 
ature  profile  at  t  —  10  s  and  t  =  600  s  due  to  transient  high-temperature  loadings.  Figure 
2.3  presents  effects  of  the  concrete  diffusivity  coefficient  a  on  transient  temperature  values  at 
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Table  2.1:  Weights  and  Abscissae  Used  in  Equation  (2.7),  N  =  1511 


k 

wk 

Xk 

1 

0.055443366310234 

0.02168694746755 

2 

0.124027738987730 

0.11268422034777 

3 

0.175290943892075 

0.27049267142189 

4 

0.191488340747342 

0.48690237038193 

5 

0.163473797144070 

0.75304368307297 

6 

0.105937637278492 

1.06093100362236 

7 

0.050027021153453 

1.40425495820363 

8 

0.016442969005267 

1.77864637941183 

9 

0.003573204214283 

2.18170813144494 

10 

4.82896509305201e-04 

2.61306084533352 

11 

3.74908650266318e-05 

3.07461811380851 

12 

1.49368411589636e-06 

3.57140815113714 

13 

2.55270496934465e-08 

4.11373608977209 

14 

1.34217679136316e-10 

4.72351306243148 

15 

9.56227446736465e-14 

5.46048893578335 

Table  2.2:  Weights  and  Abscissae  Used  in  Equation  (2.8),  M  =  1011 


k 

wk 

Xk 

1 

0.032531969510180 

0.012737849971374 

2 

0.072483896403744 

0.065802327974393 

3 

0.104004662155270 

0.156155783059660 

4 

0.121594475562980 

0.275890718366863 

5 

0.122093608318116 

0.414966322218475 

6 

0.107195747923389 

0.562009142193357 

7 

0.083077989029486 

0.704832804690269 

8 

0.056928598840185 

0.830893869740303 

9 

0.033398291993499 

0.928057569743495 

10 

0.013514893075575 

0.985992766817013 
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z  —  1  mm  and  z  —  20  mm  at  different  times.  These  results  are  consistent  with  the  graphical 
solutions  presented  by  Ju  and  Zhang.1 


Depth  (mm) 


Figure  2.2:  Concrete  Pavement  Temperature  Profile  at  t  —  10  s  and  t  =  600  s  Due 
to  Fast  Transient  Thermal  Loadings 


2.1.2  Section  Summary 

In  this  section,  rapidly  varying  1-D  temperature  profiles  in  a  homogeneous  half-space  sub¬ 
jected  to  transient  thermal  loadings  are  investigated.  The  well-known  general  solution  for 
this  problem  is  numerically  evaluated  using  efficient  Gaussian-type  integration  formulas  de¬ 
veloped  by  Steen  et  al.11  Numerical  calculations  based  on  a  uniform  initial  pavement  tem¬ 
perature  profile  and  a  model  surface  temperature  history  are  carried  out,  and  match  well 
with  published  results. 
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Figure  2.3:  Transient  Temperature  Values  at  z  =  1  mm  and  z  —  20  mm  at  Different 
Times  Due  to  Fast  Transient  Thermal  Loadings 


2.2  1-D  Temperature  Field  in  Two- layered  Pavement  Systems  Subjected  to 

High-Temperature  Transient  Loadings 

The  1-D  time-dependent  temperature  profile  in  a  homogeneous  half-space  can  be  extended 
to  a  two-layered  pavement  system,  as  shown  in  Figure  2.4.  This  idealized  two-layered  system 
can  be  eventually  used  to  analyze  a  heat-resistant  concrete  layer  over  a  conventional  concrete 
layer.  Again,  the  temperature  profile  in  a  two-layered  system  can  be  modeled  as  an  initial- 
boundary- value  problem,  where  h\  =  thickness  of  Portland  cement  concrete  (m);  h2  = 
thickness  of  the  base  layer  (m);  A j  =  thermal  conductivity  of  the  jth  layer  (kcal/m  h  °C); 
aij  =  thermal  diffusivity  of  the  jth  layer  (m2/h);  and  Tj(z,t )  =  temperature  function  for 
layer  j  (°C).  The  material  in  each  layer  is  assumed  to  be  continuous,  homogeneous,  and 
isotropic.  The  temperature  T2(z,t )  is  assumed  to  be  constant  for  z  >  H2  and  t  >  0. 
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H  i 


H2 


Figure  2.4:  Two-layered  Pavement  System 
2.2.1  Specified  Pavement  Surface  Temperature 

Similar  to  Section  2.1,  suppose  that  the  measured  transient  surface  temperature  data  is 
available,  then  1-D  time-dependent  temperature  profile  in  a  two-layered  pavement  system 
subjected  to  this  high  temperature  transient  loadings  can  be  modeled  as  the  following  initial 
boundary  value  problem 


9Tj ,  , 

Tj{z,  0)  = 
Ti(0,t)  = 
T2(H2,t)  = 

Ti(Hi,  t)  = 


=  a.; 


d2Tj 
dz 2 


%t) 


0  <  t  <  oo,  Hj_ i  <  z  <  Hj ,  j  —  1,2 


f)T 


Gj(z),  j  —  1,2  (initial  condition) 

F(t)  (first  boundary  condition) 

constant  (second  boundary  condition) 

(first  interface  condition) 

=  (second  interface  condition) 


T2{Hut) 

9 T2 


(2.11) 


dz 


where  H0  =  0,  Hi  =  hi  and  H2  —  hi  +  h2. 
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The  main  mathematical  tool  used  for  solving  the  system  (2.11)  is  a  LT.  To  facilitate  the 
derivation  of  the  solution,  we  introduce  the  variable  Uj(z,t),  j  =  1,2  below 

Uj(z,t)  =  Tj(z,t)-Tj(z,  0)  (2.12) 

For  simplicity,  we  assume  that  the  initial  temperature,  Tj(z,  0),  j  =  1,  2  is  a  constant.  From 
Equation  (2.12),  the  system  (Equation  2.11)  can  be  written  as  the  following  initial-boundary- 
value  problem: 


dui( 

at  {z’t] 

d2Ui ,  s 

=  OLi  (z,  t)  0  <  t  <  oo, 

OZz 

Ht.i  <z<Hu  i  =  1,2  (2.13) 

Ui(z,  0) 

=  0 

(2.14) 

C/r(0,t) 

=  F{t)-Tx{  0,0) 

(2.15) 

U2(H2,t) 

=  0 

(2.16) 

u^t) 

=  U2(Hut) 

(2.17) 

>■•*> 

- 

(2.18) 

and  we  assume  that  T2 ( H2 ,  t)  =  T2(H2,  0)  for  all  t  >  0. 

Let  C  denote  the  LT  operator  and  Ui(z,s )  be  the  LT  of  Ui(z,t)  with  respect  to  time  t. 
Furthermore,  the  following  operational  property12  of  LT  is  needed: 

£[f'(t)]  =  sf(s)~m  (2.19) 

where  f(s)  is  assumed  to  exist. 

Applying  LT  with  respect  to  t  to  Equation  (2.13)  in  conjunction  with  Equations  (2.14) 
and  (2.19)  yields 

muz  si  s  ^ 

]  ~  -U,(z,  s )  =  0,  Hj_x  <  z  <  j  =  1,  2  (2.20) 
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If  we  let  Vj(s) 


,  then  the  solution  of  Equation  (2.20)  is 


Uj{z,  s)  =  Aj(s)e~rjZ  +  j  =  1,  2 


(2.21) 


where  A,-(s),  Bj(s),  j  =  1,2  are  to  be  determined  using  LT  of  boundary  and  interface 
conditions.  Applying  LT  with  respect  to  t  to  the  boundary  and  interface  conditions  in 
Equations  (2.15)-(2.18)  yields  Equations  (2.22)-(2.25),  respectively: 


Ui(0,s) 

=  F(s)-C- 
s 

(2.22) 

U2(H2,s) 

=  0 

(2.23) 

UiiH^s) 

=  U2(Hi,s) 

(2.24) 

dz  {Hl’s) 

(2.25) 

where  the  constant  c  stands  for  Ti(0,  0). 

From  Equation  (2.21),  we  know 

[jfj. 

s)  =  -rjAj(s)e~rjZ  +  rjBj(s)eTjZ,  j  =  1,2  (2.26) 


Substituting  Equations  (2.21)  and  (2.26)  into  Equations  (2.22)-(2.25)  yields  the  following 
linear  system  in  which  Aj(s),  Bj(s),  j  —  1,2  are  unknown  variables 


On 

Ol2 

0 

0 

«2 1 

«22 

«23 

024 

«3 1 

«32 

«33 

034 

0 

0 

«43 

Cl44 

4i(s) 

C 1 

BiW 

0 

42(s) 

0 

0 

(2.27) 


where 
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an  —  1 


a  12  —  1 


a21  =  e-riHl  a22  =  eriHl  a23  =  -e~r2Hl  a24  =  -er2^ 

a3i  =  -Xrr1e~riHl  a32  =  Air4eri//l  a33  =  X2r2e~r2Hl  a34  =  -A  2r2er2Hl  (2.28) 

a43  =  e~r2H2  a  44  =  er2li2 

Ci  =  F(s )  -  * 

The  linear  system  (Equation  2.27)  can  be  easily  solved  by  using  Cramer’s  rule  to  give 


^i(s)  — 

Bi(s)  =  -°±I2 

4a(s)  =  2^Ainer2^  (2.29) 

B2(s)  =  -2^Xirie-r2H2 

where 


h  =  erihl  [Ayri  sinh (r2h2)  +  A 2r2  cosh (r2h2)\ 

I2  =  e~rihl  [-Ayri  sinh(r2h2)  +  A 2r2  cosh (r2/i2)] 

A  =  2  [Ayri  sinh (r2h2)  cosh(r4/ii)  +  X2r2  cosh (r2h2)  sinh(r4/ii)] 
h2  =  H2-Hx 


From  Equation  (2.4)  in  Section  2.1,  it  is  clear  that  Ts(t)  can  be  well  approximated  for 
large  t  by 

F(t)  =  285  +  49.5  ln(t),  t  >  1  (2.30) 

Since  the  LT  of  Equation  (2.30)  is  much  simpler  than  that  of  (2.4),  Equation  (2.30)  will  be 
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used  in  the  following  sample  calculation.  The  LT  of  ln(t)  takes  the  form13 


L[ln(t)] 


7  +  ln(s) 
s 


(2.31) 


where  7  ~  0.5772156  is  Euler’s  constant. 

In  view  of  Equation  (2.31),  the  LT  of  Equation  (2.30)  is 


F(s) 


285  -  49.5(7  +  ln(s)) 
s 


(2.32) 


Based  on  Tj(z,  0)  =  25,  j  —  1,2  (see  Equation  (2.5))  and  Equation  (2.32),  C\  can  be 
obtained  as 

Cx  =  -  [260  -  49.5(7  +  ln(s))l  (2.33) 

s 

Substituting  A3 (s) ,  (s) ,  j  —  1,2  in  Equation  (2.29)  into  Equation  (2.21)  and  using  the 

inverse  LT  yields 

-1  pv+ioo 

Ui(z,  t)  =  -  /  U1(z,s)estds,  0  <  z  <  Hx  (2.34) 

2vri  Jv_ioo 

1  pv+ioo 

U2(z,t )  =  —  /  U2{z,s)estds ,  H{  <  z  <  H2  (2.35) 

where  v  is  some  real  number  such  that  Uj(z,  s),  j  =  1,  2  converges  absolutely  along  the  line 
Re(s)  =  v,  where  Re(s)  denotes14  the  real  part  of  a  complex  number  s. 

Due  to  the  complexities  of  Uj(z ,  s),  j  =  1,  2,  the  closed-form  solutions  of  Equations  (2.34) 
and  (2.35)  are  difficult  to  derive,  so  we  seek  numerical  inversion  of  the  LT.  In  this  study,  the 
Gaussian  quadrature  formula15  for  evaluating  the  following  integral  of  functions  of  complex 
variables  is  employed 

-y  ru+io o  p  N  ^ 

5—  /  —  F(p)  dp  «  J2wiF(Pi)  (2-36) 

ZhI  J u — 200  P  j=l 

where  >  2  is  an  integer;  Wj,pj,j  =  1,2,  ...,1V  are  weights  and  abscissae,  respectively. 
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For  the  sake  of  completeness,  Pj,Wj ,  j  =  1,3, ...  ,9  are  listed  in  Table  2.3,  and  pj,  Wj  are 
equal  to  the  conjugate  of  Pj-i,  Wj-i  for  j  =  2, 4, ... ,  10,  respectively.15 


Table  2.3:  Abscissae  and  Weights  Used  in  The  10-point  Gaussian  Quadrature 
Formula 


* 

T 

3 

5 

7 

9 


_ Pi _ 

12.83767707781087  +  1.66606258416230b 
12.22613148416215  +  5.012719263676864* 
10.93430343060001  +  8.409672996003092i 
8.776434640082609  +  11.92185389830121* 
5.225453367344361  +  15.72952904563926* 


Wi 

-868.4606112670226  +  15457.42053305275* 
1551.634444257753  -  8439.832902983925* 
-858.6520055271992  +  2322.065401339348* 
186.3271916070924  -  253.3223820180114* 
-10.34901907062327  +  4.110935881231860* 


For  fixed  z  and  t,  let  st  =  p.  Then  complex  integrals  in  Equations  (2.34)  and  (2.35)  can 
be  written  in  the  form  of  the  integrals  in  Equation  (2.36)  as  follows: 

1  *"7+*0O  p 

Ui(z,t)  =  —  —Flip)  dp,  0  <  z  <  Hi  (2.37) 

Z'KX  J /y _ioo  P 

1  rj+ioo  p 

U2(z,t )  =  —  -F2(p)dp,  Hi  <  z  <  H-2  (2.38) 

27T*  J  ^y_lOQ  p 

where  iq (p)  =  Uj{z,  |)2  j  =  1,2.  Then  Equations  (2.37)  and  (2.38)  can  be  approximated 
using  Eq.  (2.36).  To  verify  the  validity  of  applying  Equation  (2.36)  to  evaluate  Equations 
(2.37)  and  (2.38),  a  sample  calculation  was  performed  using  parameters  given  in  Table  2.4. 

Table  2.4:  Geometry  and  Material  Parameters  Used  in  the  Sample  Calculation 


Parameters 

Value 

Layer  thickness  (m) 

hi 

0.4 

hi 

2.0 

Thermal  conductivity,  A  (kcal/m  h  °C) 

PCC  slab 

1.85 

Base  layer 

1.20 

Thermal  diffusivity,  a  (m2/h) 

PCC  slab 

0.00468 

Base  layer 

0.00360 
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In  the  sample  calculation,  it  is  assumed  that  Tj(z,  0)  =  25  °C,  j  =  1,2,  so  in  view  of 
Equation  (2.12),  the  final  solution  Tj(z,t),  j  —  1,2  is 

Tj(z,t)  =  Uj(z,t)  +  25  (2.39) 

Figure  2.5  plots  temperature  profiles  in  the  concrete  slab  at  time  t  =  10,  60,  180,  360  and 
600  s  using  temperature  solutions  for  a  two-layered  system;  Figure  2.6  illustrates  transient 
temperature  histories  from  t  —  1  to  1000  s  at  z  =  1,  10,  20,  50,  100  and  200  mm  measured 
from  pavement  surface. 


Figure  2.5:  Transient  Concrete  Slab  Temperature  Profile  for  a  Two-layered  Sys¬ 
tem  Subjected  to  Transient  Thermal  Loading 


2.2.2  Specified  Heat  Flux  from  Aircraft  Engine,  Q(t) 

In  this  case,  we  assume  that  if  the  heat  flux  emanating  from  aircraft  engine,  Q(t),  is  known, 
then  the  underlying  mathematical  model  to  estimate  the  1-D  temperature  field  in  a  two- 
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Figure  2.6:  Transient  Temperature  Values  Evaluated  at  Different  Depths  in  the 
Concrete  Slab  for  a  Two-layered  System  Subjected  to  Thermal  Loading 


layered  pavement  system,  as  shown  in  Figure  2.4,  is  given  by  the  following  equations: 


dTi 


( z,t )  = 


dt 

Tj(z,0) 

-A.ffCU) 


d2T ’■ 

aj~d^ 


[ z,t )  0  <  t  <  oo,  Hj_ i  <  z  <  Hj,  j  =  1,  2 

j  =  1,2  (initial  condition) 


=  B 


Q(t ) 


+  ^air(^)  “  Ti(0,t) 


g  -air 


T2(H2,t)  =  constant 


T!{Hut)  = 


r)T 


=  X 


T2(Hi,  t) 

ot2 


dz 


(first  boundary  conditionl(2.40) 
(second  boundary  condition) 

(first  interface  condition) 

(Hi,t)  (second  interface  condition) 


where  B  =  pavement  surface  convection  coefficient  (kcal/m2hr  °C);  Tqr  (t)  =  air  temperature 
(°C);  and  the  other  variables  are  dehned  in  Section  2.2.1.  Note  that  the  heat  input  from 
direct  solar  radiation  is  ignored  in  this  problem  due  to  the  rapid  transient  heating  of  the 
surface  by  the  aircraft  engines.  The  only  difference  between  Systems  (2.11)  and  (2.40)  is  the 
first  boundary  condition. 
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Similar  to  Section  2.2.1,  let  Uj(z,t )  =  Tj(z,t )  —  Tj(z,  0),  j  =  1,2  and  suppose  that 
Tj(z,  0)  is  a  constant,  then  the  first  boundary  condition  in  System  (2.40)  becomes 

-Ai^|(0 ,t)  =  B  ^  +  Tair(t)-T1(0,0)-C/1(0,t)  (2.41) 

Since  S>  Tair(f)  —  T\ (0,  0) ,  we  can  drop  Tair(f)  —  77(0,0)  for  simplicity.  Thus  System 
(2.40)  can  be  rewritten  in  terms  of  Uj(z,t),j  =  1,2  as  follows: 


dUj ,  , 

d2Uj ,  , 

=  ctj  -  „  (z,t)  0  <t<oo, 

oz1 

Hj-i  <  z  <  Hj,  j  —  1,2  (2.42) 

Uj{z,  0) 

=  0 

(2.43) 

=  Q(t)  -BU^t) 

(2.44) 

u2(h2,  t) 

=  0 

(2.45) 

=  U2(Hut) 

(2.46) 

- 

(2.47) 

where  it  is  assumed  that  T2(H2,t )  =  T2(H2,  0)  for  all  t. 

Similar  to  Section  2.2.1,  the  main  mathematical  tools  employed  to  resolve  the  system  in 
equation  (2.40)  are  again  a  LT  and  the  numerical  inversion  of  a  LT.  Referring  to  Section 
2.2.1,  Uj(z,s),  j  —  1,2  are  given  by  Equation  (2.21)  with  Aj(s),  Bj(s),  j  —  1,2  determined 
using  a  LT  of  boundary  and  interlayer  contact  conditions.  Applying  a  LT  to  Equation  (2.44) 
with  respect  to  t  gives 

dU 

-Ai^(0 ,s)  =  Q(s )  -  BU^s)  (2.48) 

In  the  following  sample  calculation,  Q(t),  the  step  function  representing  the  heat  flux 
emanated  from  the  aircraft  engines,  is  assumed  to  be  given  by 


I  Qo  if  h  <  t  <  t2 
Q{t)  =  l 

0  if  0  <  t  <  t\  or  t  >  t2 


(2.49) 
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where  Q 0  is  a  constant  heat  flux;  ti,  t2  are  two  time  values.  Thus,  the  LT  of  Equation  (2.49) 

is 

Q(s)  =  —  (e~stl  -  e~st 2)  (2.50) 

s 

As  in  Section  2.2.1,  A,(s),  Bj(s),  j  =  1,2  can  be  determined  using  the  linear  system 
(2.27)  with  all  the  symbols  defined  in  Equation  (2.28)  except  the  following: 


On  —  B  +  Air  i 

Oi2  =  B  -  Am 

Ci  =  —  (e~stl  -  e~st2) 

s  v  7 


and  Aj(s),  Bj(s),  j  =  1,2  are  given  as  follows: 


Ai(s) 

Bi(s) 

A2(s) 

B2(S) 


¥ 


2C^X1r1ehl(r2~ri) 

A 

C 

2  1  ^c-‘2r2h2-h1(r\+r2) 

A  1  1 


where 


(2.51) 


(2.52) 


h  =  Am  +  A 2r2  +  (A2r2  -  Am)  e  27-2,12 

I2  =  (A2r2  -  Ain)  e-27-1'71  +  (Ain  +  A2r2)  e-2^1+7-2/l2)  (2.53) 

A  =  {B  +  Am)  [Ain  +  A2r2  +  (A2r2  -  Am^-27-2'12] 

+  (B  -  Airi)e-27'lhl  [Am  -  A2r2  -  (Am  +  A2r2)e-2r2/l2] 
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Table  2.5:  Additional  Parameters  Assumed  in  The  Sample  Calculation 


Parameters 

Value 

B  (kcal/m2  h  °C) 

16.29 

Qo  (kcal/m2  h) 

90,000 

Tj(z,  0),  j  —  1,2  (°C) 

Time  variables  used  in  Q(t )  (s) 

25 

h 

10 

^2 

130 

Furthermore,  inserting  Equation  (2.52)  into  Equation  (2.21)  yields  Uj(z,s),  j  —  1,2  below 


Ui(z,s) 

=  S-  { (Am  +  A 2r2)e~riz  [l  -  e-2(rihl+r2h2~riz)] 

+  (A 2r2  -  Am)  [e~^h2+r^  -  e-2nfti+n*]  } 

(2.54) 

=  (^-2X1r1e~r^z~hl)~rihl  [1  -  e-2ra(tf2-*)j 

A 

(2.55) 

where  C\  is  given  in  Equation  (2.51). 

As  in  Section  2.2.1,  Uj(z,t),  j  =  1,  2  can  be  determined  by  an  inverse  LT  as  in  Equations 
(2.34)  and  (2.35).  The  numerical  inversion  can  be  estimated  using  the  10-point  Gaussian 
quadrature  formula  shown  in  Equation  (2.36).  In  the  sample  calculation,  the  parameters 
from  Table  2.4  are  selected  in  addition  to  those  given  in  Table  2.5. 

Figure  2.7  plots  temperature  profiles  in  the  concrete  slab  at  time  t  =  15,  30,  60,  90, 
120,  150,  and  210  s  using  temperature  solutions  for  a  two-layered  pavement  system  in  this 
section.  Figure  2.8  illustrates  transient  temperature  histories  from  t  —  1  to  1000  s  at  z  —  0, 
1,  3,  5,  8,  10,  20,  30,  and  40  mm  measured  from  pavement  surface. 

2.2.3  Sensitivity  Study 

In  this  subsection,  we  conduct  a  brief  sensitivity  study  of  the  effects  of  thermal  properties 
and  the  thickness  of  the  first  layer  on  the  temperature  profile  in  a  two-layered  system.  This 
will  give  some  clues  to  the  selection  of  appropriate  materials  having  heat-resistant  properties 
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Figure  2.7:  Transient  Concrete  Slab  Temperature  Profile  in  The  First  Layer  for  a 
Two-layered  Pavement  System  Subjected  to  Specified  Heat  Flux  from  Aircraft 


Figure  2.8:  Transient  Temperature  Values  Evaluated  at  Different  Depths  in  The 
First  Layer  for  a  Two-layered  Pavement  System  Subjected  to  Specified  Heat 
Flux  from  Aircraft 


22 

Distribution  A:  Approved  for  public  release;  distribution  unlimited. 


Table  2.6:  Geometry  and  Material  Parameters  Used  in  The  Sensitivity  Study 


Parameters 

Value 

Layer  thickness  (mm) 

hi 

60,100 

hr2 

400 

Thermal  conductivity,  A  (keal/mm  s  °C) 

Geopolymer  paste 

2.0xl0-7 

PCC  slab 

5.1xl0-7 

Thermal  diffusivity,  a  (mm2/s) 

Geopolymer 

0.2 

PCC  slab 

1.3 

for  the  surface  layer  of  airfield  concrete  pavements.  Geopolymer  materials  have  desirable 
properties  for  serving  as  an  alternative  binder  to  traditional  Portland  cement  in  producing 
paving  concrete.  These  properties  include  lower  thermal  conductivity  and  diffusivity  values, 
high  compressive  strength  at  early  age,  non-flammability,  and  high  thermal  stability.  Thus 
it  is  possible  to  construct  paving  concrete  made  from  a  geopolymer  binder  on  top  of  the 
ordinary  concrete  slab  to  limit  temperature  penetration  into  the  ordinary  concrete  layer.16 
The  following  sensitivity  study  gives  an  example  of  such  a  two-layered  system. 

The  parameters  used  in  the  sensitivity  study  are  given  in  Tables  2.5  and  2.6.  Figures  2.9 
and  2.10  plot  temperature  profiles  at  different  times,  and  Figures  2.11  and  2.12  plot  transient 
temperature  values  evaluated  at  different  depths,  for  a  two-layered  system  with  h\  =  60  mm 
and  h\  =  100  mm.  Actual  calculations  show  that  there  are  no  differences  in  the  first  nine 
significant  digits  between  calculated  temperature  values  in  generating  Figures  2.9-2.12,  i.e., 
fixing  all  the  other  parameters  and  replacing  h\  =  60  mm  by  h\  =  100  mm  does  not  change 
temperature  profiles  in  the  two-layered  system  under  the  rapidly  imposed  thermal  loading 
case.  However,  Figures  2.8  and  2.11  demonstrate  that  the  peak  temperature  values  in  the 
two-layered  system  containing  geopolymer  materials  are  lower  than  those  in  the  ordinary 
concrete  two-layered  system  at  all  depths  except  the  surface,  as  expected.  In  particular,  at 
z  =  40  mm,  the  peak  temperature  drops  from  around  100  °C  in  Figure  2.8  to  about  37  °C 
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Depth  (mm) 


Figure  2.9:  Transient  Temperature  Profile  for  a  Geopolymer-concrete  System 
(hi  =  60  mm)  Subjected  to  Specified  Heat  Flux  from  Aircraft  Operation 


Depth  (mm) 


Figure  2.10:  Transient  Temperature  Profile  for  a  Geopolymer-concrete  System 
(hi  =  100  mm)  Subjected  to  Specified  Heat  Flux  from  Aircraft  Operation 
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Figure  2.11:  Transient  Temperature  Values  Evaluated  at  Different  Depths  in  a 
Geopolymer-concrete  System  (hi  =  60  mm)  Subjected  to  Specified  Heat  Flux 
from  Aircraft  Operation 


Figure  2.12:  Transient  Temperature  Values  Evaluated  at  Different  Depths  in  a 
Geopolymer-concrete  System  (hi  =  100  mm)  Subjected  to  Specified  Heat  Flux 
from  Aircraft  Operation 
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in  Figure  2.11. 


2.2.4  Section  Summary 

In  this  section,  1-D  rapidly  varying  temperature  profiles  in  two-layered  pavement  systems 
subjected  to  transient  thermal  loadings  were  studied.  The  underlying  solution  techniques 
were  LT  and  numerical  inverse  LT.  Analytical  solutions  were  derived  for  both  the  specified 
surface  temperature  history  and  the  heat  flux  from  aircraft  engine  conditions.  Numerical 
calculations  were  carried  out  to  illustrate  the  derived  solutions.  Also,  a  brief  sensitivity 
study  of  the  effects  of  material  thermal  properties  and  the  thickness  of  the  first  layer  on  the 
temperature  profile  in  a  two-layered  system  was  conducted. 
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Chapter  3 


2-D  AXISYMMETRIC  FIELD  IN  HOMOGENEOUS  HALF-SPACE 

In  this  chapter,  analytical  solutions  of  a  2-D  axisymmetric  transient  temperature  held  are 

derived  under  the  assumption  that  the  thermal  loadings  and  surface  boundary  conditions  are 
axisymmetric.  To  take  advantage  of  axisymmetry,  a  cylindrical  coordinate  system  is  used, 
as  shown  in  Figure  3.1,  and  a  =  thermal  diffusivity  (m2/h)  and  T(r,  z,  t )  =  the  temperature 
function.  Here,  we  assume  that  the  surface  temperatures  are  available  during  the  period  of 
interest.  The  mathematical  formulation  of  this  problem  is  given  as 


dT 

~dt 

(d2T 

-  “  (  Sr*  + 

1  dT  d2T\ 

O  +02’  0  <  t  <  oo, 

r  or  ozz  J 

0  <  z  <  oo  (3.1) 

T(r ,  0,  t) 

=  F(r,t), 

(boundary  condition) 

(3.2) 

T(r,z,  0) 

=  G(r,z), 

(initial  condition) 

(3.3) 

where  F  and  G  are  assumed  to  be  continuous. 

Let  the  time  period  of  interest  be  [0,  te],  and  m  a  positive  integer.  Suppose  that  [0,fe] 
is  divided  into  2 m  sub-intervals  of  equal  length,  and  that  the  surface  temperature  at  r  =  0 
is  measured  at  two  end  points  of  each  sub-interval  except  at  time  te.  Then  the  interpola¬ 
tor  trigonometric  polynomials,  based  on  the  discrete  least  squares  approximation,  can  be 
obtained  to  approximate  F(Q,t)  as  follows17 

CL  &  m~  1  t 

F(  0,  t)  —  -+  +  cos  (mi)  +  [dk  cos  (kt)  +  bk  sin(/ct)] ,  0  <  t  <  te  —  (3.4) 

Zj  Zj  /-.III/ 

k= 1 


27 

Distribution  A:  Approved  for  public  release;  distribution  unlimited. 


with 


Figure  3.1:  Cylindrical  Coordinate  System 


t 

bk 


^  2m—  1 

—  V'  Ti  cos 
m 


1  2m—  1 

—  ^  T/  sin 
m 


'k-n  ' 

—  U  —  m) 
m 

'k-n  ' 

—  (/  —  m) 
m 


for  each  k  —  0, 1, . . . ,  m 
for  each  k  —  1, 2, . . . ,  m  —  1 


(3.5) 

(3.6) 

(3.7) 


where  XJ  =  measured  surface  temperature  at  r  =  0  at  Ith  partitioning  point  of  [0,fe],  i.e.,  at 
time  ti  =  gf- te  for  each  /  =  0, 1,  2, . . . ,  2m  —  1.  In  the  following,  F(r,  t )  is  assumed  to  have 
the  form 

F(r,  t)  =  e-^rF(0,  t)  (3.8) 


where  /i  is  some  parameter. 

For  simplicity,  we  assume  that  T(r,z,  0)  is  independent  of  z.  Also,  the  compatibilities  of 
initial  and  boundary  conditions  at  z  —  0,  t  —  0  impose  that  F(r,  0)  =  G(r,  0),  which  yields 
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T(r,  z,  0)  =  F(r,  0). 


In  view  of  the  above  discussion,  the  following  focuses  on  the  derivation  of  the  analytical 
solution  for  T(r,z,t)  satisfying  the  partial  differential  equation  (PDE)  given  in  Equation 
(3.1)  with  the  boundary  and  initial  conditions  given  by 


T(r,0,t)  =  + 


2rrm 

(1  \1 

- 1  +  7T 

—  —  nr 

4 

\2  )\ 

m—  1 

-£ 

k=  1  L 


2k7l  (\  \\  .  (2k7l 

dk  sm  I  - 1  +  7T  -  —  k  )  )  +  Ok  sill  (  — — t  —  k IT 


tp 


tp 


,  (3.9) 


where  0  <  t  <  te  (l  —  AAj ;  and 


r(r,Z,0)  =  e-'»q|  +  ^Siii 


7T  |  -  -  m 


m—  1  r 


+  £ 


k=  1  L 


ak  sin  |  7T  I  --  k 


+  bk  sin  (— kn) 


(3.10) 


where  Equation  (3.10)  is  obtained  by  setting  t  —  0  in  Equation  (3.9).  It  is  noted  that 
Equation  (3.1)  is  linear,  so  the  principle  of  linear  superposition  implies  that  the  final  solution 
satisfying  the  Equations  (3.1),  (3.9)  and  (3.10)  can  be  obtained  by  summing  up  each  solution 
satifying  Eq.  (3.1)  and  the  following  boundary  and  initial  conditions 


T(r,  0,f)  =  e  ^irA  sin(u;t  +  0)  (3.11) 

T(r,z,  0)  =  e~^rAsin(j)  (3.12) 


where  we  note  that  e~^rA  sin(u;f  +  0)  and  e~^rAsm(j)  resemble  the  variable  terms  in  the  right 
hand  side  of  Equations  (3.9)  and  (3.10),  respectively.  Hence,  the  model  initial-boundary- 
value  problem  consisting  of  Equations  (3.1),  (3.11)  and  (3.12)  will  be  considered. 

3.1  Separation  of  Variables 

The  method  of  SV  has  been  employed  to  predict  time-dependent  temperature  profiles  in 
multilayered  pavement  systems  using  the  measured  air  temperature,  solar  radiation  intensity 
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and  material  parameters.7  To  facilitate  the  derivation  of  an  analytical  solution,  the  complex¬ 
valued  function  of  real  variables,  Y{r ,  z,  t )  is  introduced,  which  is  the  solution  of  the  following 
boundary  value  problem 


dY 

fd2Y  1  DY  d2Y\ 

0  <  t  <  oo, 

0  <  z  <  oo  (3.13) 

~dt 

\  dr2  '  r  dr  '  dz2  ) 

Y(r,0,t) 

(3.14) 

Y(r ,  z,  t) 

is  bounded 

(3.15) 

where  i  is  the  imaginary  unit  number  with  i2  =  —1. 

It  is  clear  that  the  imaginary  part  of  Y(r,z,t)  satisfies  the  Equations  (3.1)  and  (3.11), 
and  in  general  does  not  satisfy  the  Equation  (3.12).  However,  the  influence  of  initial  data 
T(r,z,  0)  on  transient  temperature  distributions  gradually  decays  as  time  increases,9  and 
thus  the  solution  based  on  the  method  of  SV  can  still  give  a  reasonable  approximation  to 
temperature  at  the  point  Q(r,z,t )  for  small  z  and  large  t. 

The  following  outlines  the  main  steps  involved  in  solving  Equations  (3.13)  and  (3.15) 
based  on  the  method  of  separation  of  variables: 

1.  We  assume  that 

Y(r,  z,  t )  =  u(r,  z)ej^t+ ®  (3.16) 


then,  it  follows  that 


d  Y 
~dt 


juY 


2.  Inserting  Equations  (3.16)  and  (3.17)  into  Equation  (3.13)  yields 


jum  =  a 


d2u  Id  u  d2u\ 
dr 2  r  dr  dz 2  / 


(3.17) 


(3.18) 
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3.  We  assume  that  ur  is  bounded  at  r  =  0,  u(r,  z)  is  0(r~k )  1  and  ur(r,  z )  is  0(r~k+1 )  as 
r  — >•  oo  for  each  z  >  0  with  k  >  |.  Then  the  HT  of  order  zero  of  u(r,  2)  with  respect 
to  r,  u(£,z)  defined  below  exists19 

/»oo 

u(£,z)=  /  ru(r,z)J0(£r)dr  (3.19) 

Jo 

where  Jo(£r)  is  a  Bessel  function  of  the  first  kind,  of  order  zero. 

Applying  the  HT  on  r  to  both  sides  of  Equation  (3.18),  we  obtain  formally  the  following 
equation 

^(f,*)  ^  (£2+  W?)  u(£,z)  =  0  (3.20) 

Note  that  the  following  fact  was  used  in  deriving  Equation  (3.20)  when  u  and  ur  satisfy 
the  above-mentioned  conditions19 

r°°  f  rft  1  a  \ 

l  r  +  -rlrr)  u(r,z)Mtr)dr  =  -?u&z)  (3.21) 

4.  Solving  Equation  (3.20),  we  find  that 


u 


(f ,  z)  =  Ce~^M+jN)  +  De^M+jN) 


(3.22) 


where  M  =  yA±l;  N  =  ,  V  =  \j  1  +  ( ^2  )  ;  and  C,  D  are  constants  of 

integration  that  are  determined  using  the  boundary  condition. 


5.  The  boundedness  of  Y(r,z,t)  implies  D  =  0  in  Equation  (3.22),  and  the  inverse  HT 


1The  order  symbol  O  is  defined  as  [18,  pp.  570-571] 

f(k)  =  0[G(K)\ ,  k  — >  a  (here  a  may  be  ±00)  if 
F(k) 

absolute  value  of  ,  (  approaches  to  A  as  k  — >  a,  where  A  is  a  nonzero  constant 
G{k) 
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of  u  gives 


19 


Y(r,z,t)=  /  iC{e>e^Mz+^Nz+4,)Uir)di 
Jo 


(3.23) 


6.  Setting  z  —  0  in  Equation  (3.23)  and  considering  Equation  (3.14),  we  find  that 


19 


C(f)  =  /  rAe-»rJ0(£r)dr 

Jo 

=  A/i(£2  +  /i2)-§ 


(3.24) 


7.  Substituting  Equation  (3.24)  into  Equation  (3.23)  yields  the  complete  expression  for 
Y(r,z,t),  whose  imaginary  part,  T(r,z,t),  is  the  desired  solution  satisfying  Equations 
(3.1)  and  (3.11) 


T(r,  z,  t )  = 


f0  + 


e  Z Mz  sin  (cut  —  £Nz  +  (j))J0(^r)d^ 


(3.25) 


It  can  be  shown  that  for  fixed  t  the  improper  integral  in  Equation  (3.25)  converges 
uniformly  with  respect  to  r  and  z,  where  r  G  [  0,  oo  )  and  z  G  [  0,  oo  )  by  using  the  Weierstrass 
criterion  on  uniform  convergence  of  improper  integrals  involving  parameters.19 

In  practice,  the  improper  integral  can  be  approximated  using  numerical  integration 
schemes  such  as  Gaussian  quadrature  formulas.  However,  the  assumption  in  Equation  (3.16) 
may  not  be  valid  even  for  moderate  values  of  z  under  the  condition  of  rapidly  changing  ther¬ 
mal  loading,  for  example  z  —  40  mm  as  illustrated  in  Figure  (3.6)  below.  Therefore,  we 
propose  another  solution  method  based  on  integral  transforms  such  as  LT  and  HT. 


3.2  Integral  Transforms 

In  this  section,  we  seek  an  analytical  solution  satisfying  Equations  (3.1),  (3.11)  and  (3.12) 
based  on  LT  and  HT.  The  main  steps  involved  in  the  derivation  of  solution  are  summarized 
as  follows: 
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1.  Refer  to  Section  (2.2.1),  Let  T(r,z,s )  denote  the  LT  of  T(r,z,t )  with  respect  to  time 
t.  Applying  the  LT  with  respect  to  t  to  both  sides  of  Equation  (3.1)  yields 


s  -  A 

— T(r,z,s ) - e_Mrsin0 

a  a 


d2T  1  df  d2f 

dr 2  +  r  dr  dz 2 


(3.26) 


2.  We  assume  that  Tr  is  bounded  at  r  =  0,  T  is  0(r~k )  and  Tr  is  0(r~k+1 )  as  r  — >■  oo 
for  each  z  >  0,  t  >  0  with  k  >  |.  Applying  the  HT  of  order  zero  on  r  to  both  sides  of 
Equation  (3.26)  produces  the  ordinary  differential  equation 

-  (V  +  -)  T  =  --/u(^2  +  ii2yi  sin0  (3.27) 

dzz  V  a/  a 

where  T(£,  z,  s )  denote  the  HT  of  order  zero  of  T(r,  z,  s ). 

3.  The  solution  of  Equation  (3.27)  is 

T(£,  z,  s )  =  C(£,  s)e~^  +  D(£,  s)e> 9z  +  Tp  (3.28) 


where  /3  =  sjy  +  C(£,s)  and  D(£,s)  are  constants  of  integration,  and  Tp  stands 

for  a  particular  solution  of  Equation  (3.27)  and  is  given  by 


T  = 


Aji{fj2  +  jJ2)  3/2 
a£2  +  s 


sin  (j) 


(3.29) 


4.  Boundedness  of  T(r,z,t)  implies  that  for  fixed  complex  number  s  with  Re(s)  >  0, 
T(r,  z,  s)  is  bounded  for  r  >  0,  z  >  0,  and  it  follows  that  D(£,  s )  =  0.  Thus 


f(£,z,s)  =  C(Z,s)e-fiz  +  fp 


(3.30) 


where  C(£,s)  is  to  be  determined  using  the  BC. 
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5.  Applying  the  LT  on  t  to  both  sides  of  Equation  (3.11)  yields  T(r,0,s) 


T(r,  0,  s)  =  Ae  ,jr  (  „  U — -  cos  0  H — — -  sin  < 


Ld2  +  S 2 


ca2  +  s2 


(3.31) 


Applying  the  HT  of  order  zero  on  r  to  both  sides  of  Equation  (3.31)  gives 


T(£101s)  =  A^e  +  V2) 


2  4-  , ,21-3/2 


u  ,  s  .  , 

cos  0  H - — — -  sin  0 


cu2  +  s2 


ca2  +  s2 


(3.32) 


6.  Setting  z  —  0  in  Equation  (3.30)  and  comparing  with  Equation  (3.32)  gives 


s )  —  A/x(£2  +  n2)  3//2 


Ld 


Ld2  +  S2 


COS  0  + 


ca2  +  s2  ct£2  +  s 


sin  ( 


(3.33) 


7.  Substituting  Equation  (3.33)  into  Equation  (3.30)  and  performing  the  inverse  HT  of 
order  zero  of  T(£,  z,  s )  yields 


T(w)  =  Afi  £(£2  +  /*2r3/2 


w  ,,  1  s 

COS  0  + 


ta2  +  s2 


ta2  +  s2  cd;2  +  s 


sin  0 


ct£2  +  s 


sin  0  ^  Jo(£r)d£ 


(3.34) 


8.  Referring  to  Section  2.2.1,  the  final  solution  of  T(r,z,t )  can  be  approximated  using 
Equation  (3.34)  by  numerical  inverse  LT  methods  such  as  Gaussian-Quadrature-type 
formulas. 


3.3  Numerical  Results 

In  this  section,  numerical  results  based  on  the  above  mentioned  two  solution  methods  are 
presented.  In  the  following  calculation,  a  =  0.0035  m2/h  and  /i  =  0.01  1/m,  and  tempera¬ 
tures  at  r  —  0  with  z  —  0,  1,  5,  10,  20,  40,  60  mm  are  calculated  starting  from  t  —  0  until 
t  =  1475  s  with  an  increment  of  25  s. 
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The  surface  temperatures  at  r  —  0  for  every  25  s  from  t  —  0  until  t  =  1475  s  are  gener¬ 
ated  using  Figure  3.1  from  Jn  and  Zhang.2  Prescribed  surface  temperature  values  T(0,0,  t) 
at  t  —  0,  25, 50, . . . ,  1475  s  and  those  generated  using  the  above  mentioned  interpolatory 
trigonometric  polynomials  are  presented  in  Figure  3.2.  When  the  approximation  based  on 
the  method  of  SV  is  used,  the  surface  temperatures  T(0, 0,  t)  are  assumed  to  be  T(0,  0,  0)  at 
each  time  t  =  —300,  —275,  —250, . . . ,  —50,  —25  s  in  order  to  get  a  more  accurate  solution  of 
T(r,  z,  0)  with  r  >  0  and  z  >  0. 

For  the  numerical  results  based  on  the  method  of  separation  variables,  the  composite 
16-point  Gaussian  quadrature  formula  is  employed  to  evaluate  Equation  (3.25),  replacing 
the  upper  limit  oo  by  £  =  30,  which  is  determined  using  an  error  analysis.  The  length  of 
each  subinterval  equals  0.2  in  the  composite  Gaussian  integration  scheme,  and  the  10-point 
Guassian  quadrature  formula  used  for  resolving  inverse  LT  in  Section  2.2.1  is  employed  again 
to  obtain  numerical  values  of  T(r,  z,  t ). 

Figure  3.3  shows  the  prescribed  temperature  T(0,  0,  t)  at  7  =  0,  25,  50, ... ,  1475  s  and  the 
predicted  ones  based  on  the  methods  of  separation  of  variables  and  LT,  respectively.  Figure 
3.3  indicates  that  the  surface  temperatures  at  r  —  0  were  almost  exactly  recovered  by  the 
results  based  on  the  method  of  SV,  and  well  approximated  for  t  <  450  s  by  results  based  on 
the  LT.  The  artificial  oscillation  of  temperature  for  the  large  t  exhibited  in  the  approximation 
based  on  the  LT  is  probably  caused  by  the  error  associated  with  the  numerical  inverse  LT. 

Figures  3. 4-3. 6  present  the  predicted  transient  temperature  T(0,  z,  t )  at  z  =  1,5, 10, 20, 40, and 
60  mm  using  the  methods  of  separation  of  variables  and  LT,  respectively,  identified  in  Table 
3.1.  Figures  3. 4-3. 6  reveal  that  for  z  =  1,5,10  mm,  the  method  of  separation  of  vari¬ 
ables  gives  a  reasonable  prediction  of  temperature  except  at  small  values  of  t,  whereas  the 
method  of  LT  generates  reasonable  approximation  except  for  artificial  oscillations  exhibited 
at  t  >  700  s,  which  are  suspected  to  be  caused  by  the  numerical  inverse  LT.  For  z  =  20, 40,  60 
mm,  LT  gives  better  results  than  SV  does.  The  reason  behind  this  fact  is  that  the  latter 
does  not  use  the  initial  temperature  values  and  the  assumption  made  in  Equation  (3.16) 
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Table  3.1:  Proposed  Transient  Temperature  Prediction  Using  the  Combined  So¬ 
lution  Technique 


z  (mm) 

T{0,z,t) 

0 

based  on  SV 

1 

T(0, 1,  0),  T(0, 1,  25)  based  on  LT,  the  others  based  on  SV 

5 

T(0,  5,  0),  T(0,  5,  25),  T(0,  5, 50)  based  on  LT,  the  others  based  on  SV 

10 

T(0, 10,  0),  T(0, 10,  25)  based  on  LT,  the  others  based  on  SV 

20 

based  on  LT 

40 

based  on  LT 

60 

based  on  LT 

may  not  be  valid  in  general. 

To  take  advantage  of  the  reasonable  temperature  prediction  generated  by  each  solution 
method,  a  combined  solution  technique  is  proposed  in  this  study.  For  example,  using  the 
results  presented  in  Figures  3. 3-3. 6,  we  proposed  that  the  final  approximation  for  the  tran¬ 
sient  temperature  at  z  —  0, 1,  5, 10,  20, 40,  60  mm,  as  shown  in  Figure  3.7,  be  generated  using 
the  approaches  given  in  Table  3.1. 

3.4  Section  Summary 

In  this  section,  a  2-D  axisymmetric  temperature  field  with  specified  surface  temperature 
history  in  a  homogeneous  half-space  due  to  transient  thermal  loading  is  studied.  Two  solution 
methods  are  proposed,  one  based  on  the  method  of  SV  and  HT,  and  the  other  based  on 
LT  and  HT.  Inverse  HT  and  LT  can  be  resolved  numerically.  A  combined  approach  to 
a  solution  is  proposed  using  results  based  on  these  two  methods.  Model  calculations  show 
that  the  combined  solution  approach  gives  a  reasonable  approximation  to  the  rapidly  varying 
temperature  profile. 
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Figure  3.2:  Prescribed  T(0,  0,  t)  at  t  —  0,  25,  50, . . . ,  1475  s  and  Its  Predicted  Values 
Based  on  The  Interpolatory  Trigonometric  Polynomials 


Figure  3.3:  Prescribed  and  Predicted  Surface  Temperatures  at  r  =  0  for  Different 
Times  Based  on  LT  and  SV  Methods 
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Figure  3.4:  Predicted  Temperatures  at  r  =  0,  z  —  1  mm  and  r  =  0,  z  —  5  mm  for 
Different  Times 
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Figure  3.5:  Predicted  Temperatures  at  r  =  0,  z  =  10  mm  and  r  =  0,  z  —  20  mm  for 
Different  Times 
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Figure  3.6:  Predicted  Temperatures  at  r  =  0,  z  —  40  mm  and  r  =  0,  z  =  60  mm  for 
Different  Times 
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Figure  3.7:  Predicted  Transient  Temperatures  at  Different  Depths  Using  Results 
Based  on  Two  Methods  Described  in  This  Section 
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Chapter  4 


SUMMARY 

This  research  work  mainly  concerns  analytical  solutions  for  predicting  the  rapidly  varying 
temperature  profile  in  pavements  subjected  to  transient  thermal  loading.  Analytical  1-D 
and  axisymmetric  2-D  solutions  are  developed  in  this  technical  report. 

The  general  1-D  solution  for  a  pavement  temperature  profile  in  homogeneous  half-space 
with  specified  surface  temperatures  is  well-known  and  can  be  resolved  numerically  using 
some  powerful  Gaussian-quadrature  type  integration  formulas  recommended  in  this  report. 
Analytical  solutions  for  the  rapidly  varying  temperature  profile  in  two-layered  pavement  sys¬ 
tems  are  systematically  investigated  in  this  study  and  it  is  shown  that  they  can  be  used  to 
analyze  the  thermal  effect  of  an  innovative  heat-resistant  concrete  layer  overlying  a  conven¬ 
tional  concrete  base  layer.  The  main  mathematical  tools  employed  in  deriving  temperature 
profiles  in  two-layered  pavement  systems  are  LT  and  numerical  inversion  of  LT.  One  dimen¬ 
sional  solutions  are  derived  for  input  conditions  both  of  specified  surface  temperature  and  of 
heat  flux  from  aircraft  engines.  Model  calculations  suggest  that  the  derived  1-D  analytical 
solutions  can  capture  the  rapidly  changing  transient  pavement  temperature  profile  in  both 
homogenous  half-space  and  two-layered  systems. 

For  the  2-D  axisymmetric  Dirichlet  problem  (i.e.,  pavement  surface  temperatures  are 
known)  in  homogeneous  half-space,  specified  axisymmetric  transient  surface  temperatures 
are  assumed  and  two  solution  methods  are  developed.  The  first  method  is  based  on  HT  and 
the  method  of  SV,  and  the  other  is  based  on  LT  and  HT.  Numerical  experiments  suggest  that 
a  solution  combining  results  generated  by  these  two  analytical  methods  can  give  reasonable 
predictions  for  actual  rapidly  varying  temperature  profiles. 
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Appendix 


This  appendix  lists  the  main  Matlab  and  Fortran  codes  developed  in  this  study. 

1.  1-D  Temperature  Field  in  a  Homogeneous  Half-Space  (Matlab  codes) 

•  homo _varied _z.m:  calculate  temperature  profile  at  a  fixed  time  t 

•  homo_varied_t.m:  calculate  temperature  history  at  a  particular  depth 

2.  1-D  Temperature  Field  in  a  Two-layered  Pavement  System  (Matlab  codes) 

(a)  Specified  Pavement  Surface  Temperatures 

•  twolayer_varied_z.m:  calculate  temperature  profile  in  the  concrete  layer  at  a 
fixed  time  t 

•  twolayer_varied_t.m:  calculate  temperature  history  at  a  particular  depth  in 
the  concrete  layer 

(b)  Mixed  Boundary  Conditions,  i.e.,  Specified  Heat  Flux  from  Aircraft  Engine 

•  twolayer  mixed  be  varied  z.m:  calculate  temperature  prohle  in  the  concrete 
layer  at  a  fixed  time  t 

•  twolayer_mixed_bc_varied_t.m:  calculate  temperature  history  at  a  particular 
depth  in  the  concrete  layer 

3.  2-D  Axisymmetric  Temperature  Field  in  a  Homogeneous  Half-Space  Subjected  to  Spec¬ 
ified  Pavement  Surface  Temperatures  (Fortran  source  codes) 

(a)  Separation  of  Variables  Method 

•  temphomo2.f:  main  code  to  calculate  temperature  history  at  a  particular 
location  in  the  concrete  layer 

•  gauss. f  and  wl.f:  two  subroutines 
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(b)  Laplace  Integral  Transformation  Method 


•  homo  It.  f:  main  code  to  calculate  temperature  history  at  a  particular  location 
in  the  concrete  layer 

•  gauss. f  and  wl.f:  two  subroutines 
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LIST  OF  SYMBOLS  AND  ABBREVIATIONS 


1- D 

2- D 
B 
BC 
F 
G 

Hi 

H2 

HT 

MZr) 

LT 

Q(t) 

sv 

T 

Tj 

f(Z,z,s) 

Uj 

U3 

Y(r,z,  t) 

hi 

h2 

3,M,N 

r 

s 


one- dimensional 
two-dimensional 

pavement  surface  convection  coefficient  (kcal/m2hr  °C) 
boundary  condition(s) 

function  appeared  in  the  boundary  condition 
function  appeared  in  the  initial  condition 
hi 

hi  +  h2 

Hankel  transform 

first  kind  of  Bessel  function  of  order  zero 
Laplace  transform 

heat  flux  emanating  from  the  aircraft  engines  (kcal/nr2  h) 
separation  of  variables 
temperature  function  (Celsius  degree) 
temperature  function  for  layer  j  (Celsius  degree) 

Laplace  transformation  of  T(r,  z,t)  with  respect  to  t  followed  by  Hankel  transfor¬ 
mation  on  r 

shifted  temperature  function  for  layer  j  (Celsius  degree) 

Laplace  transformation  of  Uj 

complex- valued  function  whose  imaginary  part  is  the  desired  temperature  function 
thickness  of  Portland  cement  concrete  (m) 
thickness  of  the  base  layer  (nr) 
integers 

radial  variable  (nr) 

Laplace  transformation  variable 
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t  temporal  variable  (s) 

Wk  weights  used  in  the  Gaussian  quadrature  type  formula 
Xk  abscissae  used  in  the  Gaussian  quadrature  type  formula 
2  depths  measured  from  pavement  surface  (m) 

a  thermal  diffusivity  of  material  (m2/h) 

a,j  thermal  diffusivity  of  the  j th  layer  (m2/h) 

A j  thermal  conductivity  of  the  jth  layer  (kcal/m  h  °C) 

H  parameter  to  account  for  the  temperature  variation  along  the  radial  direction 

(l/m) 

£  Hankel  transformation  variable 

(j)  phase  angle  (non-dimensional) 
uj  frequency  (1/hour) 
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